%% Fig12b_Replication.m
% 复刻论文 Fig.12(b)：RM-DCSK 在固定 SNR 下的 BER 与 Spread Factor N 的关系
% 调试步骤：
% 1. 每种映射生成混沌序列后进行归一化，并绘制前 1000 点波形，检查其统计特性
% 2. 检查 RM-DCSK 编码后信号的均方功率是否一致
% 3. 固定 SNR（如 25 dB），对不同 N 进行仿真，统计 BER，并用 semilogy 绘制 BER vs. N
%
% 注意：请确保参数（如映射参数、帧结构、归一化处理）与论文描述一致。

clear; clc; close all;
rng(0);  % 固定随机种子，确保结果可重复

%% 1. 参数设置
SNR_dB_fixed = 25;         % 固定 SNR（dB）
N_values = 10:10:100;        % Spread Factor N 取值范围（例如 10、20、…、100）
numBits = 1e4;             % 发送比特数；调试时可先用1e5，调试通过后可提高至5e5以上

% 定义待测试的混沌映射名称（2D 映射 + 1D 映射对比）
chaosMapNames = { '2D-MLM','2D-MSM','2D-MTM','2D-SMM',...
                  'Logistic','Sine','Tent' };

% 各映射对应参数，确保与论文一致
paramsMap = { struct('mu',0.2,  'k',2.08, 'x0',0.2, 'q0',0.1), ... % 2D-MLM
              struct('mu',-0.1, 'k',1.88, 'x0',0.8, 'q0',0.6), ...  % 2D-MSM
              struct('mu',-0.6, 'k',1.78, 'x0',0.8, 'q0',-0.8), ...  % 2D-MTM
              struct('mu',0.5,  'k',2.33, 'x0',1,   'q0',1), ...    % 2D-SMM
              struct('mu',3.68, 'x0',0.2), ...                      % Logistic (1D)
              struct('mu',0.93, 'x0',0.2), ...                      % Sine (1D)
              struct('mu',1.49, 'x0',0.2) };                        % Tent (1D)

% 帧结构：每个比特对应 2*N 个采样点，
% 前 N 点为参考段，后 N 点为数据段（数据段乘以比特值）

%% 2. 主仿真循环：对不同 Spread Factor N 计算 BER
numMaps = length(chaosMapNames);
BER_matrix = zeros(numMaps, length(N_values));

for mapIdx = 1:numMaps
    mapName = chaosMapNames{mapIdx};
    params = paramsMap{mapIdx};
    fprintf('Simulating map: %s\n', mapName);
    
    for iN = 1:length(N_values)
        N = N_values(iN);
        
        % 生成混沌序列（总长度 = 2*N*numBits）
        chaosSeq = generateChaosSequence(mapName, numBits, N, params);
        
        % Debug：绘制混沌序列前 1000 点波形
        if iN == 1 && mapIdx == 1  % 仅对第一映射的第一个 N 绘制一次
            figure('Name',['Chaos Sequence: ', mapName]);
            plot(chaosSeq(1:min(1000, end)));
            title(['前 1000 个采样: ', mapName, ', N=', num2str(N)]);
            xlabel('采样点'); ylabel('幅值');
            fprintf('Map: %s, N=%d, 归一化前均值=%.4f, 标准差=%.4f\n', mapName, N, mean(chaosSeq), std(chaosSeq));
        end
        
        % 对混沌序列进行归一化：使均值为0, 标准差为1
        chaosSeq = (chaosSeq - mean(chaosSeq)) / std(chaosSeq);
        
        % Debug: 输出归一化后均值、标准差
        if iN == 1 && mapIdx == 1
            fprintf('After normalization: mean=%.4f, std=%.4f\n', mean(chaosSeq), std(chaosSeq));
        end
        
        % 生成随机比特（-1, +1）
        dataBits = randi([0,1], 1, numBits)*2 - 1;
        
        % RM-DCSK 编码
        txSignal = RM_DCSK_Encode(dataBits, chaosSeq, N);
        
        % Debug: 检查 txSignal 均方功率
        txPower = mean(txSignal.^2);
        if iN == 1 && mapIdx == 1
            fprintf('Map: %s, N=%d, txSignal 均方功率=%.4e\n', mapName, N, txPower);
        end
        
        % 通过 AWGN 信道（固定 SNR）
        rxSignal = awgn(txSignal, SNR_dB_fixed, 'measured');
        
        % RM-DCSK 解码
        rxBits = RM_DCSK_Decode(rxSignal, N);
        
        % 计算 BER
        BER = sum(rxBits ~= dataBits) / numBits;
        BER_matrix(mapIdx, iN) = BER;
        
        % Debug: 输出部分 SNR 下的 BER 信息
        if iN == 1  % 对每个映射，N=第一个值时输出信息
            fprintf('Map: %s, SNR=%d dB, N=%d, BER=%.4e\n', mapName, SNR_dB_fixed, N, BER);
        end
    end
end

%% 3. 绘制 BER vs. Spread Factor N 图 (半对数坐标)
figure('Name','Fig12(b) Replication - BER vs. Spread Factor');
hold on; grid on;
markers = {'o','s','^','d','+','x','*'};
lineStyles = {'-','-','-','-','-','-','-'};

for mapIdx = 1:numMaps
    semilogy(N_values, BER_matrix(mapIdx,:), lineStyles{mapIdx}, ...
        'Marker', markers{mapIdx}, 'LineWidth',1.5, 'MarkerSize',8, ...
        'DisplayName', chaosMapNames{mapIdx});
end
xlabel('Spread Factor N');
ylabel('BER');
title(['固定信噪比', num2str(SNR_dB_fixed), ' dB','下 误码率随扩频因子 N 的变化']);
legend('Location','Best');

%% ============== 子函数 ==================

%% 产生混沌序列函数
function chaosSeq = generateChaosSequence(mapName, numBits, N, params)
    % 每个比特帧长度 = 2*N, 总长度 = 2*N*numBits
    totalLen = 2 * N * numBits;
    chaosSeq = zeros(1, totalLen);
    
    switch mapName
        case '2D-MLM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DMLM(x, q, mu, k);
            end
        case '2D-MSM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DMSM(x, q, mu, k);
            end
        case '2D-MTM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DMTM(x, q, mu, k);
            end
        case '2D-SMM'
            mu = params.mu; k = params.k; x = params.x0; q = params.q0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                [x, q] = iteration_2DSMM(x, q, mu, k);
            end
        case 'Logistic'
            mu = params.mu; x = params.x0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                x = mu * x * (1 - x);
            end
        case 'Sine'
            mu = params.mu; x = params.x0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                x = mu * sin(2*pi*x);
            end
        case 'Tent'
            mu = params.mu; x = params.x0;
            for i = 1:totalLen
                chaosSeq(i) = x;
                if x < 0.5
                    x = mu * x;
                else
                    x = mu * (1 - x);
                end
            end
        otherwise
            error('Unknown chaos map name!');
    end
end

%% RM-DCSK 编码函数
function txSignal = RM_DCSK_Encode(dataBits, chaosSeq, N)
    % 每比特对应一帧，帧长 = 2*N
    L = length(dataBits);
    txSignal = zeros(1, 2*N*L);
    idx = 1;
    for i = 1:L
        ref = chaosSeq(idx : idx+N-1);                % 参考段
        dat = dataBits(i) * chaosSeq(idx+N : idx+2*N-1);% 数据段，乘以比特
        txSignal((2*N*(i-1)+1):(2*N*i)) = [ref, dat];
        idx = idx + 2*N;
    end
end

%% RM-DCSK 解码函数
function rxBits = RM_DCSK_Decode(rxSignal, N)
    % 每帧 2*N 点，采用相关器判决
    L = length(rxSignal) / (2*N);
    rxBits = zeros(1, L);
    idx = 1;
    for i = 1:L
        ref = rxSignal(idx : idx+N-1);
        dat = rxSignal(idx+N : idx+2*N-1);
        correl = sum(ref .* dat);
        if correl > 0
            rxBits(i) = +1;
        else
            rxBits(i) = -1;
        end
        idx = idx + 2*N;
    end
end

%% 迭代方程：2D-MLM
function [x_next, q_next] = iteration_2DMLM(x, q, mu, k)
    x_next = mu * x * (1 - x) + k * cos(q) * x;
    q_next = q + x;
end

%% 迭代方程：2D-MSM
function [x_next, q_next] = iteration_2DMSM(x, q, mu, k)
    x_next = mu * sin(2*pi*x) + k * cos(q) * x;
    q_next = q + x;
end

%% 迭代方程：2D-MTM
function [x_next, q_next] = iteration_2DMTM(x, q, mu, k)
    if x < 0.5
        x_next = mu * x + k * cos(q) * x;
    else
        x_next = mu * (1 - x) + k * cos(q) * x;
    end
    q_next = q + x;
end

%% 迭代方程：2D-SMM
function [x_next, q_next] = iteration_2DSMM(x, q, mu, k)
    x_next = mu * x + k * cos(q) * x;
    q_next = q + x;
end
